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Abstract. The dynamics of a small Prandtl number binary mixture in a laterally heated cavity is studied 
numerically. By combining temporal integration, steady state solving and linear stability analysis of the 
full PDE equations, we have been able to locate and characterize a codimension-three degenerate Takens- 
Bogdanov point whose unfolding describes the dynamics of the system for a certain range of Rayleigh 
numbers and separation ratios near S=-l. 



1. Introduction 

Double-diffusive fluxes occur when convection is driven by thermal and concentration gradients, and the 
temperature and con centration diff usivities take different values. This phenomenon has relevance for nu- 
merous applications (Turner, 1985), and from a theoretical point of vie w presents very interesting dynamics, 
including chaos JCross and Hohenber 2 119931 iKnobloch et all Il986h . We are interested here in the case 
of horizontal gradients dTurnerlll980Hjiang et all Il99lh . In~this configuration, quiescent (conductive) so- 
lutions can exist when thermal and solutal buoyancy forces exactly compensate each other. This occurs 
only for a very particular value of the separation ratio (S — — 1, see below), but has allowed for theoretical 
analysis by studying the stability of the conductive s olution JGhoraveb and MojtabiUl99"7tlXin et allll997t 
iBardan et all l200(A iBergeon and Knoblochl l2002h . In recent work, we address ed this case for a small 
Prandtl number binary mixture, including only the Soret effect (Mec a~e7aiT. 120 04). Results showed a quite 
interesting bifurcation scenario by varying the Rayleigh number. In particular, we found an orbit that is 
born in a global sadd le-loop bifurcatio n, becomes chaotic in a period doubling cascade, and disappears in a 
blue sky catastrophe dShilnikovlfl997ll . This orbit is the only stable solution in a large interval of Rayleigh 
numbers. In this paper we analyze this system in greater depth to determine how the scenario associated to 
the origin of this orbit is modified when the value of S is changed to larger (i.e. less negative) values. It is 
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relevant to assess to what extent the dynamics depends on tuned values of the parameters, or whether it is 
fairly robust to these changes. Moreover, by varying a second parameter we gain access to a richer portrait 
of the system, obtaining bifurcation lines and points of codimension two. 

We have numerically integrated the full PDE equations in a region near S = — 1, combining steady state 
solving, numerical continuation, linear stability analysis, and temporal integration. The results show that 
distinct bifurcations of the S = — 1 case (namely two saddle nodes and a global saddle loop) approach each 
other in the region near S = —0.9 until only a Hopf bifurcation remains, in a scenario consistent with the 
unfolding of a codimension-three degenerate Takens-Bogdanov point. 

The outline of this paper is as follows. In Section|2]we detail the model and the numerical procedure. In 
Section|3]the behavior of the system for S — —1 is reviewed. In Section|4]we extend these results by letting 
both Ra and S vary. Finally, the discussion of the results and some concluding remarks are presented in 
Section |5] 



2. Basic equations and numerical methods 

We consider the 2-D flux of a binary mixture in a rectangular cavity Q of length d and height h. The 
aspect ratio T = d/h has been chosen to be 2. The cavity is laterally heated, maintaining different constant 
temperatures at the opposed vertical boundaries. AT is the difference between both temperatures. On the 
horizontal boundaries, a linear temperature profile is imposed. All the boundaries are taken to be no-slip and 
with no mass flux. In these conditions the dimensionless equations in Boussinesq approximation explicitly 
read 

d t u + (u • V)u = -VP + crV 2 u 

+ aRa[{\ + S) (-0.5 + x/T) + 9 + SC]z, 
d t + (u ■ V)0 = -v x /T + V 2 6>, (1) 
d t C + (u ■ V)C = -v x /T - tV 2 {6 - C), 
V ■ u = 0. 

In these equations lengths, times and temperatures are scaled with h, t K = h 2 /n and AT, respectively, k 
being the thermal diffusivity. u = (v x ,v z ) is the (dimensionless) velocity field in (x, z) coordinates, P is 
the pressure over the density, 9 and C are deviations from a linear horizontal profile of the temperature and 
of the rescaled concentration of the heavier component, respectively. The dimensionless parameters are the 
Prandtl number a = v/n, the Rayleigh number Ra = agh? AT /vk and the Lewis number t = D/k, where 
v denotes the kinematic viscosity, g the gravity level, a the thermal expansion coefficient, and D is the 
mass diffusivity. The separation ratio is defined by S = C (l — Co)- St, where St is the Soret coefficient, 
Co is the actual value of the concentration of the heavier component in the homogeneous mixture, and [3 
is the mass expansion coefficient (positive for the heavier component). Finally, boundary conditions are 
written as 

u = = n- V(C-0) = 0, at df2. (2) 

Thes e boundary conditions are no t compatible with t h e transformation of the Soret equations into those use d 
bv lGhoraveb and Moitabil < 19971) : Ixin et alJ < 19971) : iBardan et aD (EoOO): Ber geon and Knoblochl J2002I) . 
Note also that Eqs. Q together with boundary conditions (J2ji, are invariant under a rotation tt around the 
point (r/2, 1/2 ) as (x, z) — > (T — x, 1 — z), (v x ,v z ,9, C) — > (—v x , —v z , —9, — C). Therefore the system is 
Z 2 -equivariant ( Kuznet sovl Il998l) . From now on solutions invariant (non-invariant) under n will be called 
symmetric (non-symmetric). 

Eqs. and boundary con ditions have been solved by a second order time-splitting algorithm, 
proposed by Hug ues and Randriamampianinal Jl998l> . applied to a pseudo-spectral Chebyshev method. 
To calculate steady solutions, we have adapted a pseudoespectral first-o r der ti m e-stepping formulatio n 
to carry out Newton's m ethod, as described by Ma mun and Tuckermanl l ll995l) : iBergeon et alJ 1^998); 
Ixin and Le Ouerd J200 lb . In the preconditioned version of New ton's iteration , the c orresponding linear 
system is solved by an iterative technique using a GMRES package jFravsse et all E003). The linear stability 
analysis of the steady states is conducted by computing the leading eigenvalues of the Jacobian by means 
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Figure 1. Stationary solutions diagram for S = —1 varying _Ra, as obtained by Mec a et all 120041) . Continuous lines: stable states, 
dashed lines: unstable states. 



of Arnoldi's method, using routines from the ARPACK package. For numerical calculations the chosen 
parameters have been Prandtl number a = 0.00715 and Lewis number r = 0.03. The system has been 
discretized in space by using 72 x 48 and 90 x 60 mesh points in steady calculations, giving both resolutions 
equivalent results. For example, increasing the resolution, the Rayleigh number of the turning points varies 
less than a 0.1% Temporal integration was used basically to follow orbits with very long periods, and in 
particular in regimes with divergent periods where a fit of the divergence itself was needed. Thus, much 
more computation time was required. A mesh grid of 60 x 30 points has proved to be sufficient to obtain 
results that did agree fairly well with steady calculations in the cases where the comparison applied. 



3. Scenario for S = — 1 

In this section we summarize the behavior of the system for a separation ratio S = — 1 JMeca etaill2004l) . 
The bifurcations diagram of the steady solutions is shown in Fig.^ In this figure, the Nusselt number Nu, 
defined as the ratio of the heat flux through the hot wall to that of the corresponding conductive solution, is 
represented as a function of the Rayleigh number Ra. For the sake of clarity, only one of the non-symmetric 
solutions related by tt is represented. For this value of S the conductive solution is allowed. For small 
Ra this solution is stable, but loses stability at Ra — 541.9 through a transcritical bifurcation (T\). The 
resulting solutions are symmetric by the rotation n, and are characterized by a central main roll accompanied 
by secondary ones in opposed corners. The supercritical branch is stable only until Ra — 542.4, where a 
small non-symmetric branch connects it to the c onductive solutions through pitchfork bifurcations, a case 
analogous to that reported by Bardan et alJd2000l) . We center the discussion here on the solutions originating 
from the subcritical branch. This solution is stabilized by a saddle node bifurcation at Ra = 99 (SN{) 
and loses stability again in a Pitchfork bifurcation at Ra = 245 (P), where a stable non-symmetric branch 
appears. In Fig. |2 (left) we see that the breaking of the symmetry confines the main roll to one lateral 
side. Continuing the symmetric branch, a supercritical Hopf bifurcation at Ra = 2137 (Hi) (maintaining 
symmetry) is found. At Ra = 2253 the bifurcating periodic solution gains stability in a Pitchfork bifurcation. 
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Figure 2. Non symmetric steady solution for S = — 1. Left: Stream lines for the stable solution at Ra = 888. Right: Detail of the 
non-symmetric branch in the region of the two turning points, indicating the sign of the real part of the leading eigenvalues. The open 
circle denotes a global saddle-loop bifurcation. 



Furthermore, the steady non-symmetrical solution undergoes a saddle-node bifurcation at Ra — 1209 
(SN 2 ) and has a saddle-saddle turning point at Ra = 865.6 (TP). Further along this branch, a Hopf 
bifurcation can be found at Ra = 2218 (i?2)- A detail of this branch in the region of the two turning points 
is represented in Fig. |2 (right). Here, the variable Ek, related with the kinetic energy and defined as 
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Ek = -= / / v-vdxdz, (3) 

1 Jx=0 J z=0 

is plotted versus the Rayleigh number Ra. In this figure we have also included the sign of the real part of 
the leading eigenvalues; a parenthesis is used to indicate a complex conjugated pair. Between both turning 
points an homoclinic saddle loop connection at Ra = 1184 (SL) gives birth to a periodic non-symmetric 
solution ( Me ca et ail l2004|) . This orbit, which is the only stable solution in a wide range of values of the 
Rayleigh number Ra, is characterized by very long periods and a spiking behavior. In particular at the SL 
connection its period diverges logarithmically as expected. 

When the Ray leigh number is increased, a very interesting complex behavior of this orbit arises 
jMeca et ail I2004I) . Firstly, at Ra = 2137 the orbit starts to show ripples, reflecting the frequency corre- 
sponding to the Hopf bifurcation Hi, while the period increases dramatically. In the region near Ra = 2235 
the orbit undergoes a period doubling cascade, becoming chaotic. At Ra = 2257.5 the chaotic a t tractor 
disappears in a blue sky catastrophe, in a scenario similar to that proposed by S hilnikov and Turaevl (2000) 
in which both length and period of an orbit diverge at the bifurcation point. 



4. Results for S > -1 

We have performed both steady state and temporally dependent calculations of the system for different 
values of S above —1. We have centered our research on the solutions from which the non-symmetric orbit 
is born in a global saddle loop (SL) connection. Namely, we are referring to the non-symmetric branch 
with two successive turning points (SN2 and TP), between which the SL connection is found. This is the 
situation represented in Fig. |2] (right) for S = — 1. 

This configuration of bifurcations changes when S is increased. The sequence of events in this region is 
quite complex, but the final situation is simple. Remarkably, for values around S — —0.8920 only a Hopf 
bifurcation remains. To analyze this process we present in Fig.|3]results for this branch and different S values. 
In this figure, branches and local bifurcations are found by continuing steady solutions and performing linear 
stability analysis, whereas global bifurcations are located by temporal integration. To do this, we fix a value 
of S, and starting from a periodic solution as initial condition, we vary slightly the Rayleigh number while 
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Figure 3. Detail of the bifurcations diagram for several values of S > — 1. Kinetic energy is represented versus Ra. Full circles denote 
Hopf bifurcations. Global bifurcations (either SL or SNIC) are represented by open circles. 



monitoring the value of the period of the final stable orbit. Then, the connection is located at the point 
where the period diverges. 

The first qualitative change occurs between S = — 1 and —0.9, where the global connection SL moves 
toward the saddle node SN 2 , becoming a SNIC (saddle node on an invariant circle) in a saddle node 
loop (SNL) codimension-two bifurcation. This change is manifested in the law of divergence of the period 
T found by temporal integration of the bifurcating orbit for Ra toward the homoclinic connection. The 
divergence for the SL is logarithmic, i.e. 

T~ -- - log (Ra-Ra c )+ A, (4) 
A 

where Ra c is the value for which the global connection takes place and A is the eigenvalue of the jacobian 
matrix associated to the unstable direction of the hyperbolic solution. The divergence changes to square 
root for the SNIC: 

T ~ , B + A. (5) 
V Ra — Ra c 

Here, Ra c corresponds to the position of the saddle node. Furthermore, as we increase S, the point at 
which the two positive eigenvalues merge to become a pair of conjugated complex values (see Fig. fright) 
approaches the turning point (TP) until a Takens-Bogdanov (TB) codimension-two bifurcation occurs there 
for S near —0.8990. In Fig|3]we see the branch for S = —0.897, for which the TB has unfolded to a Hopf 
bifurcation and a global saddle loop connection, which constitute the birth and the end of an orbit that has 
been calculated by temporal integration. The divergence of the period of this orbit at the saddle loop is 
again logarithmic, as in Eq. |4j Therefore for this value of S we obtain both divergences, logarithmic for 
the SL and square root for the SNIC. We have performed fits of the periods for both points. Indeed, the 
A value obtained by the fit, A/it = 0.0173, agrees fairly well with the eigenvalue obtained by the stability 
calculation of the steady state (A = 0.0166). For the SNIC, Ra c {fit) = 856.64, and the value of SN 2 is 
Ra c = 856.67. In this situation the eigenvalue that became positive at SN 2 gains stability at TP. At the 
Hopf bifurcation the solution loses stability again. 

When S is slightly increased, an additional codimension-two SNL bifurcation is found at SN 2 , the 
SNIC becoming an SL moving away from the saddle node. At this moment, both SL are approaching each 
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Figure 4. Periods of the homoclinic orbits found at S = —0.8966 near their corresponding SL, together with the corresponding 
logarithmic fits. Left: divergence located at Ra = 854.853; Right: divergence located at Ra = 855.595 



other, a situation represented in the branch S — —0.8966 in Fig. [3] The divergences of the periods of the 
corresponding homoclinic orbits are shown in Fig. 0] They are both logarithmic, with fitted values for the 
eigenvalues being Xf it = 0.00968 and 0.0118, according well with the steady calculations A = 0.00820 and 
0.0134, respectively. 

Very soon afterwards they touch each other and disappear. By S = —0.894 no global connection remains 
in this branch. It is now possible for the two turning points (SN2 and TP) to annihilate each other in a 
codimension-two cusp bifurcation. That occurs for S = —0.8928. The final situation, in which only the 
Hopf bifurcation is found, is represented by the S = —0.892 branch of Fig.|5J 



5. Discussion and concluding remarks 

In the preceding section we have numerically studied the change of an interesting branch of non-symmetric 
solutions when a second parameter (the separation ratio S) is varied. The scenario that emerges from these 
results consists of a series of codimension-two bifurcations arranged to enable two destabilizing saddle 
nodes of a branch connected to an orbit by a global bifurcation to disappear, resulting in a simpler situation 
with only a local Hopf bifurcation. 

These results can be summarized by drawing bifurcation lines and codimension-two points in an (S, Ra) 
plot. This is shown in Fig. |5] In this figure the two turning point lines meet in a cusp bifurcation at 
S = —0.8928. We can track the larger eigenvalues of the system along these lines. We find that at S = —0.8990 
the solution on the line TP has a double zero eigenvalue. This indicates the presence of a Takens-Bogdanov 
(TB) codimension-two point. Indeed, this can also be checked by noting that, by performing linear stability 
analysis, a branch of Hopf bifurcations, with very low frequency starting from zero, also begins at this 
point. It should be stressed that in order for both turning points to be effectively annihilated TP should 
become stabilizing. That happens at the TB point. 

The existence of a TB bifurcation also implies the existence of a branch of homoclinic connections 
(saddle-loop). This is, in fact, one of the few analytical methods of proving the existence of a homoclinic 
orbit. We have used the time evolution code to situate the line of homoclinic connections as detailed in Sec. 
|4] and checked if the period of the orbit diverges following the logarithmic law (Eq. |4j. Indeed the values 
found for A fu did agree well with the eigenvalues calculated from the stability of the steady solutions. The 
fitting of the value of Ra c yields the points marked in Fig.|5] We see that the line that joins these points is 
connected with the saddle-node line corresponding to the SN 2 points. This connection is SNLi (located 
very close to S — —0.897), one of the saddle-node separatrix-loop points l lSchecteru l987) mentioned above. 
It is at this point that the branch of homoclinic connections merges with the saddle-node line to become 
an SNIC line. This is confirmed by the change of the logarithmic law for the divergence of the period (Eq. 
|4j to the square root of Eq. |5j 
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Figure 5. Bifurcation lines in the (S, Ra) parameter space (see text). 



If we keep following the SNIC line toward more negative separation ratios, we find the other codimension- 
two saddle-node separatrix loop point (SNL 2 , at S = —0.92, not shown in Fig.|5j where the SNIC line 
ends, giving birth again to a homoclinic bifurcations line separated from the saddle-node line. Following 
this line until S=-l we recover the SL point at Ra = 1184. 

This particular configuration of codimension-two bifurcations is very far from being peculiar of this 
problem. It is found in many other areas, such as bursting oscillat ions in neural o r biological system s 
)de Vriesifl99^:lBorisuklfl99l . population dynamics (Bazvkin, 1985), laser dynamics (Mav ol et all l2QQ2). 
and some diffusively coupled systems jKanamami and SekineL 120031) . Theoretically, all this behavior can be 
reproduced as a 2-dime nsional slice in paramete r space of the unfolding of a codimension-three degenerate 
Takens-Bogdanov point foumortier et allfl99ll) . also known as the Dumo rtier-Roussarie-So t omay or (DRS) 
bifurcation. In particular, we are referring to the focus case described bv lDumortier et all dl99ll) . This is 
also consiste nt with the fact that the presence of the SNIC cannot be related to a local codimension-2 
phenomenon (Golubistk v~ alll200lT) . 

Th e scenario in which the DRS bifurcation appears can be described by a planar vector field ( D umortier et all 
119911) . i.e. it is essentially two-dimensional. This suggests the possibility of calculating the normal form 
coefficients of the bifurcation from the PDE's, which is an interesting problem from a theoretical point of 
view. 
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